"""
PDEJulia.jl - A comprehensive Julia package for solving partial differential equations

This package provides a complete ecosystem for solving PDEs in Julia, inspired by the 
py-pde Python package but optimized for Julia's performance and features.

# Main components:
- Grids: Spatial discretization (Cartesian, spherical, cylindrical)
- Fields: Data structures on grids (scalar, vector, tensor)
- PDEs: Pre-defined and custom PDE implementations
- Solvers: Time integration methods (explicit, implicit, adaptive)
- Storage: Data persistence and analysis
- Trackers: Real-time monitoring and analysis during simulation
- Visualization: Plotting and animation tools
- Tools: Utilities for configuration, parameters, and helpers

# Quick start:
```julia
using PDEJulia

# Create a 2D Cartesian grid
grid = CartesianGrid([[0, 10]], [64])

# Define initial condition
u0 = ScalarField(x -> sin(π * x[1] / 10), grid)

# Define diffusion PDE
pde = DiffusionPDE(diffusivity=0.1)

# Solve for 10 time units
result = solve(pde, u0, t_span=10.0)
```

# Export structure:
The package exports all necessary functions and types for typical PDE workflows.
Advanced users can access internal modules directly.
"""
module PDEJulia

# Core Julia packages
using LinearAlgebra
using SparseArrays
using StaticArrays
using Logging

# Differential equations
using DifferentialEquations
using OrdinaryDiffEq
using BoundaryValueDiffEq

# Symbolic and modeling
using Symbolics
using ModelingToolkit
using ForwardDiff
using ChainRulesCore

# Data handling
using HDF5
using JLD2
using CSV
using DataFrames

# Numerical tools
using FFTW
using Interpolations

# Visualization
using Plots
using PlotlyJS
using Makie

# Utilities
using ProgressMeter
using BenchmarkTools
using PrecompileTools

# Version information
const VERSION = v"1.0.0"

# Include all submodules
include("grids/grids.jl")
include("fields/fields.jl")
include("pdes/pdes.jl")
include("solvers/solvers.jl")
include("storage/storage.jl")
include("trackers/trackers.jl")
include("visualization/visualization.jl")
include("tools/tools.jl")

# Import modules
using .Grids
using .Fields
using .PDEs
using .Solvers
using .Storage
using .Trackers
using .Visualization
using .Tools

# Core API - Grid types
export 
    # Abstract types
    AbstractGrid, AbstractCartesianGrid, AbstractSphericalGrid,
    
    # Concrete grid types
    UnitGrid, CartesianGrid, PolarSymGrid, SphericalSymGrid, CylindricalSymGrid,
    
    # Grid functions
    grid_laplace, grid_gradient, grid_divergence, grid_curl,
    get_boundary_conditions, set_boundary_conditions

# Core API - Field types
export 
    # Abstract types
    AbstractField, AbstractDataField,
    
    # Concrete field types
    ScalarField, VectorField, Tensor2Field, FieldCollection,
    
    # Field operations
    gradient, laplacian, divergence, curl, integral, norm,
    interpolate, apply_operator

# Core API - PDE types
export 
    # Abstract types
    AbstractPDE,
    
    # Concrete PDE types
    PDE, DiffusionPDE, WavePDE, AllenCahnPDE, CahnHilliardPDE,
    KPZInterfacePDE, KuramotoSivashinskyPDE, SwiftHohenbergPDE,
    
    # PDE functions
    evolution_rate, solve_steady_state,
    solve_laplace_equation, solve_poisson_equation

# Core API - Solver types
export 
    # Abstract types
    AbstractSolver,
    
    # Concrete solver types
    ExplicitEulerSolver, ImplicitEulerSolver, RungeKuttaSolver,
    AdamsBashforthSolver, CrankNicolsonSolver, AdaptiveSolver,
    
    # Solver functions
    solve, solve!, step!, Controller,
    registered_solvers

# Core API - Storage
export 
    # Storage types
    MemoryStorage, FileStorage, HDFStorage,
    
    # Storage functions
    save_field, load_field, get_time_series, extract_time_series

# Core API - Trackers
export 
    # Tracker types
    ProgressTracker, PlotTracker, SteadyStateTracker,
    DataTracker, CallbackTracker,
    
    # Tracker functions
    track!, get_tracker_data

# Core API - Visualization
export 
    # Plot functions
    plot, plot!, animate, animate!,
    plot_field, plot_time_series,
    
    # Movie functions
    create_movie, save_animation

# Core API - Tools
export 
    # Configuration
    Config, get_config, set_config!,
    
    # Parameters
    Parameter, ParameterList,
    
    # Utilities
    parse_duration, get_package_version,
    spherical_to_cartesian, cartesian_to_spherical

# Advanced workflow functions
export 
    # High-level workflows
    run_simulation, parameter_sweep, stability_analysis,
    
    # Auto-configuration
    auto_solver, auto_tracker, recommend_grid,
    
    # Integration helpers
    create_workflow, setup_simulation

"""
    solve(pde::AbstractPDE, initial_state, t_span; solver=nothing, kwargs...)

High-level interface for solving PDEs.

# Arguments
- `pde::AbstractPDE`: The PDE to solve
- `initial_state`: Initial condition (field or array)
- `t_span`: Time span (number or tuple)
- `solver`: Solver to use (auto-selected if nothing)
- `kwargs...`: Additional arguments passed to solver

# Returns
- Solution trajectory or final state depending on storage options

# Examples
```julia
# Simple diffusion
grid = UnitGrid([64])
u0 = ScalarField(x -> exp(-(x[1]-0.5)^2/0.1), grid)
pde = DiffusionPDE(diffusivity=0.1)
result = solve(pde, u0, 10.0)

# With custom solver and tracking
solver = AdaptiveSolver(rtol=1e-6)
tracker = PlotTracker(interval=1.0)
result = solve(pde, u0, (0.0, 10.0); solver=solver, tracker=tracker)
```
"""
function solve(pde::AbstractPDE, initial_state, t_span; 
               solver=nothing, tracker=nothing, storage=nothing, kwargs...)
    # Auto-select solver if not provided
    if solver === nothing
        solver = auto_solver(pde, initial_state)
    end
    
    # Setup default storage if not provided
    if storage === nothing
        storage = MemoryStorage()
    end
    
    # Setup controller
    controller = Controller(solver; tracker=tracker, storage=storage, kwargs...)
    
    # Run simulation
    return controller(pde, initial_state, t_span)
end

"""
    run_simulation(pde, initial_state, t_span; config...)

Comprehensive simulation runner with automatic configuration.

This is a high-level interface that handles common simulation patterns
including parameter sweeps, stability analysis, and performance optimization.
"""
function run_simulation(pde::AbstractPDE, initial_state, t_span; 
                       auto_optimize=true, parallel=false, kwargs...)
    if auto_optimize
        # Optimize grid resolution, solver, and parameters
        config = optimize_simulation_config(pde, initial_state, t_span)
        return solve(pde, initial_state, t_span; config..., kwargs...)
    else
        return solve(pde, initial_state, t_span; kwargs...)
    end
end

"""
    parameter_sweep(pde_template, parameters, initial_state, t_span; kwargs...)

Run parameter sweep over multiple PDE configurations.

# Arguments
- `pde_template`: Function that creates PDE given parameter dict
- `parameters`: Dict or NamedTuple of parameter ranges
- `initial_state`: Initial condition
- `t_span`: Time span
- `kwargs...`: Additional solver arguments

# Returns
- Results dictionary with parameter combinations as keys

# Examples
```julia
# Sweep over diffusivity values
pde_template = params -> DiffusionPDE(diffusivity=params.D)
parameters = (D = [0.1, 0.5, 1.0, 2.0],)
results = parameter_sweep(pde_template, parameters, u0, 10.0)
```
"""
function parameter_sweep(pde_template, parameters, initial_state, t_span; 
                        parallel=false, kwargs...)
    # Implementation would handle parameter combinations and parallel execution
    results = Dict()
    
    # Generate parameter combinations
    param_combinations = generate_parameter_combinations(parameters)
    
    if parallel
        # Use distributed computing for parameter sweep
        results = parallel_parameter_sweep(pde_template, param_combinations, 
                                         initial_state, t_span; kwargs...)
    else
        for params in param_combinations
            pde = pde_template(params)
            result = solve(pde, initial_state, t_span; kwargs...)
            results[params] = result
        end
    end
    
    return results
end

# Helper functions for auto-configuration
function auto_solver(pde::AbstractPDE, initial_state)
    # Analyze PDE characteristics and recommend solver
    if is_stiff(pde)
        return ImplicitEulerSolver(dt=0.01)
    elseif requires_high_accuracy(pde)
        return AdaptiveSolver(alg=Tsit5(), rtol=1e-8)
    else
        return ExplicitEulerSolver(dt=0.01)
    end
end

function optimize_simulation_config(pde::AbstractPDE, initial_state, t_span)
    # Analyze problem and optimize configuration
    config = Dict()
    
    # Analyze grid resolution requirements
    if needs_fine_resolution(pde, initial_state)
        config[:adaptive_refinement] = true
    end
    
    # Analyze time stepping requirements
    if is_stiff(pde) || has_fast_dynamics(pde)
        config[:solver] = ImplicitEulerSolver()
    else
        config[:solver] = AdaptiveSolver()
    end
    
    return config
end

# Precompilation
@compile_workload begin
    # Precompile common workflows
    grid = UnitGrid([32])
    u0 = ScalarField(x -> exp(-x[1]^2), grid)
    pde = DiffusionPDE(diffusivity=0.1)
    solve(pde, u0, 1.0)
end

end # module PDEJulia